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On the identification of substructure in phase-space using 
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ABSTRACT 

We study the evolution of satellite debris to establish the most suitable space to 
identify past merger events. We confirm that the space of orbital frequencies is very 
promising in this respect. In frequency space individual streams can be easily iden- 
tified, and their separation provides a direct measurement of the time of accretion. 
We are able to show for a few idealised gravitational potentials that these features 
are preserved also in systems that have evolved strongly in time. Furthermore, this 
time evolution is imprinted in the distribution of streams in frequency space. We have 
also tested the power of the orbital frequencies in a fully self-consistent (live) N-body 
simulation of the merger between a disk galaxy and a massive satellite. Even in this 
case streams can be easily identified and the time of accretion of the satellite can be 
accurately estimated. 

Key words: galaxies: formation - galaxies: kinematics and dynamics - methods: 
analytical - methods: N-body simulations 



1 INTRODUCTION 

In the current standard cosmological model, known as 
ACDM, galaxies like our own Milky Way are formed bottom- 
up, through the merger and accretion of smaller building 
blocks th at come together du e to their gravitational attrac- 
tion (e.g. lWhite fc Reeslll978l ). 

This model is being continuously tested and many 
questions and puzzles remain to be addressed. An exam- 
ple is the large number of bound substructures predicted 
to orbit galaxies like the Milky Way and M31 compared 
to the o bserved abundance of satellite galaxie s in these 
systems (|Klvpin et al.1 1 19991 : iMoore et all 1 19991). although 
sever al solutions have been proposed (e.g. Bullock et al.l 
I2OOOI ). The formation of realistic disk galaxies in full cos- 
mological simulations rema ins another great challenge (e.g. 
ISteinmetz fc Navarrol 19991) yet s teady progress is being 
made (e.g. iGovernato et all |2004| . 120071 : IScannapieco et al] 
120091 '). 

One strong test of the current paradigm may be per- 
formed through observations of the stellar halos of galaxies 
like the Milky Way. For example, if the Galaxy was formed 
in a hierarchical fashion via mergers, then we should be able 
to find fossil signatures of these events in the present day 
phase-space distribution of halo stars. This idea was strongly 
ex alted after the disco very of the disrupting Sgr dwarf galaxy 
by llbata et akl l|l994 ). In the following years further debris 
from accretion events was discovered, not only in the halo 
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but also in the disk of our Galaxy (Eggen"'l996'; 'ibata et alj 
2003;_Yanny ct al. 2003; Navarro, Hclmi & Freeman 200i; 
[Helmi et aLll2006l ). although some of this substructure may 
be of dynamical origin (|Antoia et al.l l2009l : iMinchev et al.l 
120091 ). 

However, if the hierarchical scenario is correct, we are 
still far from having unveiled every single stream associ- 
ated to each merger event the Galaxy ha s experienced in its 
Ufetime. iHelmi. White fc Springel (|2003l ) by combining nu- 
merical simulations with analytic work predicted that there 
should be several hundreds of col d stellar streams present i n 
the vicinity of the Sun (see also IVogelsberger et al.|[2008l ). 
The short dynamical time scales, especially in the Solar 
neighbourhood, are the basic reason behind why such a large 
number is expected. An accreted satellite, orbiting in the 
inner regions of the Galaxy will give rise to multiple stellar 
streams in a re l atively short period of time. As shown by 
iHelmi fc Whit3 { 19991 ) the spatial density of these streams 
is a strongly decreasing function of time. Consequently, sub- 
structure associated to a past accretion event is expected to 
have a very low-density contrast, making its detection very 
difficult from its distribution in space. On the other hand, 
due to the conservation of phase-space density, the velocity 
dispersion of a stream will decrease as time goes by, making 
the streams tighter in the velocity domain. 

The arguments above highlight the need for full 6D 
phase-space information to completely disentangle the pre- 
dicted wealth of substructure. Furthermore, samples of 
at least 1,000 halo stars would appear to be necessary 
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for this enterprise. Statistical arguments using the revised 
NLTT proper motion survey (New Luyt en Two Tenths, 
iGould fc Saiimll2003l : ISalim fc Gouldl 12003') have been used 
to measure the amount of granulari ty (associated to moving 
groups) in the stellar halo near Sun. lGouldll|2003l ) has shown 
that no streams are present in this region of the Galaxy that 
contain more than approximately 5% of the halo stars near 
the Sun. This estimate is in good agreement with the ex- 
pectations described above, but does not provide a direct 
confirmati on of the predictions . Spectroscopic surveys such 
as R AVE (|Zwitter et al.1 120081 ) and SEGUE (|Yannv et al.1 
|2009|'1 should be helpful in identifying further nearby streams 
(e.g. Ismith et al1l2009l : iKlement et"ai1l2009l ') but the break- 
throu gh will take place onl y with the astrometric satellite 
Gaia (|Perrvman et al.ll2001^ . Gaia, expected to be launched 
in 2011, will measure the positions and motions of 10^ stars 
with very high accuracy. 

Given the datasets that soon will become available it is 
natural to ask what is the optimal way to discover merger 
debris. More precisely, the question we would like to address 
in this work is what is the best space to look for substructure 
if our goal is to disentangle how the Galaxy was assembled. 

The first steps in this direction have been presented 
bv lHelmi fc de Zeeuwl l|2000l ). These authors analysed sim- 
ulations of the disruption of satellite galaxies in a fixed 
Galactic potential and concluded that the "integrals-of- 
motion" space defined by the energy, the total angular 
momentum and its z-component {E,L,Lz) is very suit- 
able for identifying the remains of past accretion events. 
Even in fully self-consistent cosmological simulations, dis- 
rupted satellites still appear as coherent struc t ures in this 
"integ rals of motion" space l|Knebe et al.ll2005^ . iFoiit" et al.l 
l|2006l ) demonstrated the power of complementary chemical 
abund ance information to disentangle debris from different 
events. lArifvanto fc Fuchi |2006) have used for samples of 
nearby stars the space defined by {U^ + 2V^y^^ and V, 
where U is the radial velocity in the Galactic plane, and V is 
the ve locity component in the dir ection of rotation. More re- 
cently iMcMininiABinne^ (|2008l ') (hereafter MB08) showed 
that action-angles coordinates, and in particular, the orbital 
frequencies could conform a very convenient set of variables 
to identify nearby substructure. 

In this paper we explore further the power of the orbital 
frequencies. We start our journey by briefly reviewing the 
concept of action-angle variables in Section 2. In Section 

3 we discuss in depth the space of orbital frequencies. We 
characterise the distribution of debris in a few simple cases 
and in particular, in a static spherical potential. In Section 

4 we analyse how the structure of this space is altered in 
a time dependent (still rigid) potential. Finally, in Section 

5 we focus on a live N-body simulation of the accretion of 
a massive satellite onto a pre-existing thin disk. We discuss 
and summarise our results in Section 6. 



2 THE ACTION-ANGLE VARIABLES 

Let us consider a dynamical system with a time- independent 
Hamiltonian H{p, q), where (p, q) is a set of 2n canonical 
coordinates. The dynamical state of this system is governed 



by the Hamilton's equations 

. _ dH . _dH 

dq-i dpi 



(1) 



with i = 1, . . . , n. 

The evolution of a system may be expressed much more 
simply by performing a canonical transformation of coordi- 
nates such that the equations of motion are 



P, = -^ = 0, 



dH{P) 
dPi 



with solutions 



f^>(P), (2) 



(3) 



If the system under study is such that the motion is periodic, 
then this set of canonical variables is known as action-angle 
variables. 

In a spherical potential 4>(r), the radial action is 

Jr = - f ^ dr-J2[E-$(r)]r'^ -L'' (4) 
in r 

where L is the total angular momentum, ri and r2 the turn- 
ing points in the radial direction and 



Jtt> — Lz 



Je — L — Lz 



(5) 



are the angular actions. The frequencies of motion Q.i — 
dH/dJi, can be derived by diff'erentiation of the implicit 
function 
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^Jr-- / ' dri\/2[S-$(r)]r2-L2 (6) 



defined by Eq. (31). For more details, see e.g.. 
iGoldsteTnT Poole fc SaSol (|200lf ). 



3 FREQUENCY SPACE FOR STATIC 
POTENTIALS 

3.1 Toy-Model 

The space of orbital frequencies in the context of isolating 
debris from accreted satellites was first introduced by MB08. 
For completeness we here review its main characteristics. 
We also refer the reader to Section 4 of MB08. For this 
purpose we have developed a simple toy-model and followed 
the evolution of a satellite in this space. 

Let us consider a satellite (i.e. a system of particles that 
are initially strongly clustered in phase-space) orbiting un- 
der the influence of an arbitrary gravitational potential. For 
simplicity, we consider a system with only two independent 
frequencies of motion. We will also assume that the initial 
distribution in action-angle space of the satellite follows a 
multivariate Gaussian with dispersions of ag — 7r/4 and 
a J = 7r/3 in the angles and actions, respectively (and no 
correlation among the different directions). In this space, 
the location of a given particle will vary only due to the 
evolution of the angle coordinates. The rate of evolution, i.e. 
the frequencies, is completely determined by the underlying 
potential in which our satellite evolves (see Eq. (|3}). The 
other two coordinates, i.e., the actions, will remain constant 
in time. 
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Figure 1. Distribution of particles in actions (upper panels) and frequency (lower panels) space located in a small region of the angle 
space at three different times. Each clump represents a particular stream crossing this region at the time considered. The solid lines 
indicate fli = est whereas the dotted-dashed lines represent lines of constant Qi, as shown in the bottom panels. 



Here we shall assume that we can express the frequen- 
cies in the following way: 

fli = Ci Ji + C2J2, 

(7) 

= C3J1 + C4 J2, 

where d are constants. In general, the relation between the 
frequencies and the actions is more complex, and depends 
on the specific form of the gravitational potential. 

Typically, we have access to samples of stars located in 
a particular region of space such as, for example, the So- 
lar neighbourhood. This constraint acts as a filter, implying 
that we only observe the subset of stars with orbital prop- 
erties such that they are in this particular location at the 
time of observation. To mimic this, we focus on the par- 
ticles located in a small region of the angles' space in our 
toy-model. Figure [l] shows the distribution of these particles 
in the space of actions and of frequencies at three different 
times (0 < ti < t2 < ts). 

As discussed by MB08 the particles passing through a 
particular location of (angle-)space are distributed in differ- 
ent patches or streams. With time, the number of streams 
increases while their extent in frequency space decreases. 
When the particles cross our location for the first time, they 
have not had enough time to become significantly spread in 
angle. We then observe only one big stream (left panel of 
Fig. [TJ. At some later time, particles with the highest fre- 



quencies begin to overtake the slowest ones and therefore 
we observe multiple streams (since the angles are periodic) 
each with its own characteristic frequency. The size of each 
individual stream at any given time is dependent on the size 
of the small region under consideration and on the initial 
spread of the system. 

Since the orbital frequencies are constant (for a time 
independent gravitational potential) the location of a par- 
ticle in frequency space remains unchanged. The changing 
pattern shown in Fig. [T] is simply a consequence of the fact 
that we restrict our analysis to particles at a given spatial 
location. The gaps thus correspond to particles which are 
located outside our window. 

The network of patches is shaped by lines of constant 
frequency fii = est. As a consequence, and because of the 
relation implied by Eq. ([7}, the distribution of particles in 
action space is also very regular. However, in general the 
relation between frequencies and actions will not be linear, 
and therefore the distribution of particles in action space is 
typically more complex (see e.g. Figure [S] below) . 

The number of clumps observed along each frequency 
direction at a given location can be estimated by considering 
the initial spread in the frequency. For example, at fl2 = est, 
AQi = n™'"' - n™'". After some time t = is, the angle 6*1 
between the particles with the largest and smallest Qi has 
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increased to 



A6»i(t3) = Anit3 + Ae»i(to) 



(8) 



with A6li(to) the initial separation at time to. Since 6\ is a 
27r-periodic variable the number of streams A^'i at ^2 = est 
is 



iVi = 







' Aflit^ + AOiitoY 






2-K 




2n 




2n 



(9) 



— , i = l,2. 



(10) 



where the last approximation holds if ta is long enough such 
that Afiita >> A6'i(to). 

Another quantity of interest is the distance between two 
adjacent streams in frequency space, 

The time t may be considered as the time since disruption 
of our system (i.e. its constituent stars have been tidally 
stripped, and move only under the influence of the host 
gravitational field). Eq. (llOp states that this time can be di- 
rectly estimated by measuring SQi. Moreover, systems will 
have different values of Sfli depending on the time since dis- 
ruption, which implies that this characteristic scale could be 
used to recognise the streams originating in objects plausi- 
bly accreted at different epochs. 

3.2 Spherical potentials 

We now consider dynamical systems with an underlying 
spherical potential, <I>(r). Due to the conservation of angular 
momentum L, the motion of stars are constrained to remain 
in a plane, implying that there are only two independent 
frequencies of motion: Qr associated to the radial oscillation 
and associated to the angular oscillation in this plane. 

3.2.1 Characterising frequency space 

Our previous analysis has shown that the space of frequen- 
cies appears to be well suited to identify streams. Here we 
will characterise this space and how its properties depend 
on the form of the host gravitational potential. 

The angular and radial periods of an orbit are related 
to the respective frequencies through 



T(P = 27r/f^0, Tr = 27V /rir. 



(11) 



Once a radial oscillation is completed the azimuthal angle 
has increased by an amount A(/!> — fl^Tr or, equivalently 



thus 



Tr 



2-K 
A^' 



(12) 



In general 2ti/ A(j> will not be a rational number and 
therefore the orbit will not be closed. However there are 
two gravitational potentials in which all orbits are closed, 
those of an homogeneous sphere and of a point mass. No- 
tice that any other spherical system may be considered to 
have a mass distribution between these limiting cases, i.e. 
its density gradient will be neither as shallow nor as steep 
as in these two cases. Orbits in a time-independent homo- 
geneous mass distribution have A^ = vr, whereas for the 
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Figure 2. Top: Distribution of possible orbits in a Plummer 
sphere in frequency space. The solid line corresponds to Qr = 
2n^, as in the case of orbits in a homogeneous sphere. The dashed 
line shows Or- = O0, as found for all orbits in a point mass distri- 
bution. The left wedge of the distribution corresponds to circular 
orbits. Bottom: Relation between the radial frequency Or- and the 
orbital apocentre-to-pericentre ratio. 



Kepler problem A(/> = 27r. From Eq. (|12|) this implies that 
the frequencies of a bound orbit in any spherical mass dis- 
tribution are constrained to lie inside the subspace defined 
by f2r/2 ^ Vl^ fir. The extent to which this subspace is 
probed by orbits depends on the exact distribution of mass. 

To examine this, let us consider the Plummer potential, 
$, generated by a density p, where 

$(r) = 



Vr2 + 62 



p(r) 



3M 



1 + 



-5/2 



(13) 



47r63 V o^ 

IPlummejllQllI ). We sample the frequency space in this po- 
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and density p, are given by 
$(r) = 

p(r-) = 



-GM 
r + b' 

bM 



(14) 



Figure 3. Distribution of possible orbits in a Hernquist profile 
(black dots) in frequency space. For comparison, the distribution 
for a Plummer sphere is shown in grey. Solid and dashed lines 
represent as before the natural limits of this space, i.e. fir = 20^ 
and ilr = ^tt> respectively. Note that because the Hernquist inner 
density profile is cusped orbits with high frequencies do not bend 
over to the boundary associated with the homogeneous sphere. 



tential by randomly drawing 10* orbits from the distribution 
function associated to the Plummer sphere with parameters 
M = 10^^ M© and 6 « 22 kpc. The upper panel in Fig. 
shows the distribution of these orbits in frequency space. 
The solid line and dashed lines correspond to Q,. — 2Q,/, (ho- 
mogeneous sphere) and fir — (point mass), respectively. 
Note that orbits with small frequencies cover the whole re- 
gion between these two lines. Instead, orbits with high fre- 
quencies tend to bend over towards the line given by the ho- 
mogeneous sphere. This is because particles on orbits with 
low frequencies populate the outer regions of the system, 
and hence perceive a potential close to that of a point mass. 
On the other hand those on high frequency orbits are gener- 
ally confined to the central regions of the system. Since the 
Plummer density profile has a core these particles will feel a 
potential resembling that of the homogeneous distribution. 
The strict cut-off at a value of fir ~ 42 Gyr~^ therefore cor- 
responds to the characteristic frequency associated to the 
homogeneous mass distribution, i.e for r << fe, a particle in 

the Plummer model should have a period of Tr ~ \/^^ 
which in our case corresponds to 0.15 Gyr, or fir ~ 42 
Gyr-\ 

For a given fl^, (nearly) circular orbits will be those 
with the largest radial period, i.e. the smallest fir frequency. 
Thus such orbits are located on the wedge seen in the top 
panel of Fig. (2] This is also illustrated in the bottom panel of 
the same figure, where we have plotted the ratio of apocentre 
to pericentre distances of orbits as a function of fir- As we 
can see, orbits with the lowest fir for a given fl^ = est have 

'^apo/^per — 1. 

We now consider a system with a cuspy density dis- 
tribution, namely the Hernquist profile whose potential, <1>, 



27rr (r -I- 6)3 

iHernquistll 19901 ). We set M = 10^^ Mq and b ^ 22 kpc. As 
before, we sample the available phase-space with 10* orbits. 
The black points in Fig. |3] show the resulting distribution 
of orbits in frequency space. For comparison, we overplot- 
ted the frequencies for the Plummer profile (in grey). At the 
low frequencies end both distributions overlap as expected 
(point mass regime) . However there is a significant difference 
in the region populated by orbits with high frequencies. Due 
to the lack of a core in the Hernquist density distribution, 
the bend over towards fir — 2fl^, i.e. the regime of the ho- 
mogeneous mass distribution, has disappeared, and there is 
no upper limit to the radial or angular frequency of motion. 



3.2.2 Comparing different projections of phase- space 

Various spaces have been propos ed to identify debris from 
disrupted satellites. For example, iHelmi fc de Zeeuwl (|2000l ) 
introduced the space of energy, total ang ular moment u m and 
its 2-component {E,L,Lz). Later on. iHelmT et al.l (|2006l ) 
showed that substructure can be identified by looking at the 
apocentre, pericentre and space. In this space, streams 
are represented by extended structures, roughly along a line 
of constant eccentricity. In this section we will compare the 
distribution of debris in these and other previously used 
spaces to that in frequency space. To this end we consider 
the accretion of a satellite galaxy onto a spherical host. 

We represent the host with the Plummer profile dis- 
cussed in the previous section. The satellite is assumed to 
be spherical and to have an isotropic velocity ellipsoid, and 
is modelled with a multivariate Gaussian in 6D with a dis- 
persion of (Jx ~ 1 kpc and CTv = 22 km/s~^. We set the 
central particle of the satellite on an orbit with apocentre 
T'apo ~ 35 kpc, pericentre rpcr ~ 8 kpc and inclined by 56 
degrees with respect to the z-axis. We follow then the evolu- 
tion of the satellite for approximately 10 Gyr. For simplicity, 
in this experiment the self-gravity of the satellite is not con- 
sidered. As a consequence, we expect to find a larger number 
of streams at any given epoch, with a larger spread in or- 
bital parameters, in comparison to the self-gravitating case. 
This is because particles can drift apart from each other 
right from the start, whereas when their self-gravity is con- 
sidered, this will only happen during successive pericentric 
passages where they are (typically) stripped from the parent 
satellite. 

In Fig. |4] we show the spatial distribution of satellite 
particles at three different times. To study the distribution 
of particles in various projections of phase-space we focus 
on those inside an sphere of radius 4 kpc centred at r « 10 
kpc at the final time, as shown in the rightmost panel of this 
Figure. 

The results are shown in Fig. [5] for six different spaces. 
Grey points represent all the particles of the satellite 
whereas black points only those inside the sphere of interest. 
In the top left corner (panel a) we show their projection in 
frequency space. As in the toy model discussed in Section 
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Figure 4. X-Y distribution of particles from an accreted satellite at three different times, as indicated on eacli panel. The black circle 
shows the location of a "Solar neighbourhood" sphere. 



13.11 for a given fi^ = cat we find multiple streams in fir- 
Note that there are typically more streams at constant Q.^ 
than at constant 57^- This is because Tr ^ Ttf, ^ 2Tr, and 
hence mixing occurs faster in the r-direction. 

Notice that some of the groups of particles found at 
given Qr and « est may be decomposed into smaller 
structures. This can be understood from the radial velocity 
vs. radial distance plot shown in the bottom right panel 
of Fig. [5] Streams with pericentres inside the sphere under 
consideration, appear in frequency space as a single compact 
structure. For the rest we observe two separate structures, 
each of these associated to particles located just before or 
after their corresponding pericentre. 

The top right panel of Fig. [S] shows the distribution of 
particles in a projection of velocity space. Although several 
features are visible, the individual streams are clearly less 
well separated than in panel a). This is also the case for the 
radial velocity vs. radial distance plot (panel f). 

The two middle panels of Fig. [5] show the distribution 
of particles in the E — (left) and apocentre vs. pericen- 
tre (right) spaces. The number of structures in the E — 
space is smaller and slightly less sharp than in frequency 
space. On the other hand, streams are well detached from 
each other in the apocentre vs. pericentre space. Finally, in 
the bottom left panel we show the distribution of particles 
in the space of actions, Jr vs. J^. As in the E — Lz space, 
distinction of the various streams is less clear in this space, 
and hence it appears to be less suitable for finding individ- 
ual debris streams. However, it is important to note that 
we are only considering two of the three available actions. 
By projecting the particle's distribution into the principal 
directions of action space it is also possible to find a well 
defined distribution of streams (see Figure 7 of MB08). 

This analysis shows that frequency space is a very suit- 
able space to identify streams associated to accretion events. 
Moreover, as we explained in Section 13.11 the way streams 
are distributed in this space can allow us to estimate quan- 
tities such as time since disruption. This last characteristic 
makes the space of orbital frequencies more appealing than 
the other spaces previously considered. 



3.2.3 Estimating the time of accretion 

As discussed before the separation between adjacent streams 
along each direction of the space of orbital frequencies yields 
a direct estimate of the time since disruption, e.g. a satellite 
accreted 10 Gyr ago should give rise to streams separated 
by a scale 6Q,r ~ 5S^,^> = 27r/10 Gyr~^. 

By considering information from the distribution of par- 
ticles not only in frequency but also in angle space MB08 
developed a powerful method that, by a process of iteration 
over a guess potential, allows to estimate time of disruption 
and to pin down the true underlying potential. 

In this section we explore a different method which, only 
based on the regularity seen in frequency space, provides an 
accurate estimate of time since disruption. As an example 
we focus on the streams shown in panel a) of Fig. [S] Be- 
cause there is a unique scale (all streams are separated by 
the same amount) a straightforward way to compute the 
separation between adjacent streams is to use Fourier anal- 
ysis techniques. We proceed by first creating an image from 
the scatter plot in frequency space, to which we then apply 
a Fourier transform. In the final step we compute the power 
spectrum. To obtain an image from our data we grid the 
frequency space with a regular N x N mesh of bin size A in 
both ilr and directions. In this way, the number of par- 
ticles in the {nr,n^) bin of the grid is h{nr,n^). We apply 
a two dimensional discrete Fourier transform to this image: 

JV-l JV-l 

H{kr,k^) = exp{2nikrnr / N) 

nT-=0 n^— 

exp{2TTik^n^/N)h{nr,n^), (15) 

with kr, kgy = . . . , ^ - 1. 

In the top panel of Fig. [S] we show the image of the 
panel a) of Figure [5] in Fourier space. This has been com- 
puted on a grid of 200 by 200 elements. The axes represent 
the wavenumbers in each direction {kr/{NA), k^/{NA)) 
while the colour coding show the square value of each of 
the Fourier components \H{kr, ktf,)\^ . As we can see from 
this image, and as expected from our previous discussion, 
a regular pattern is present. To measure the characteristic 
scale we consider two different "slits" (of 4 bins width) of 
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400 




Figure 5. Comparison of the various spaces commonly used to identify merger debris. Grey dots show the distribution of all satellite 
particles whereas the black dots represent the particles inside the sphere shown in Fig. |4] As in the toy model example discussed in 
Section |3. II streams in frequency space (top left panel) are distributed in a regular pattern. Notice that the distribution of particles in 
frequency space is much better defined than in any other of the spaces considered. 
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the Fourier image around kr — and — 0. These regions 
are delimited by the black lines in the top panel of Fig. |6l 
We can now compute the one-dimensional power spectrum 
of the image along each direction as 



log(IH(k,k^)l^) 
5 10 15 



P(0) = —\H{0,0)\ 



Pikr) = [\H{kr,0)f + \H{-kr,0)f] 



(16) 



for fc, = l,...,(f -1), 

P{N/2)^^\H{-N/2,0)\^ 

and analogously for P{k^). Recall that the one-dimensional 
power spectrum we compute is, in practice, an average over 
the slits considered along the kr and k^ directions. 

The results are shown in the middle and the bottom 
panels of Figure |S1 for P(kr) and P{k^), respectively. The 
scales present in the power spectrum P{kr) are associated 
to separations in the fir direction (since they correspond to 
k^ — 0). The large power for wavenumbers close to zero, i.e. 
P(0) represents the average power over the whole image. All 
other peaks in the spectrum are related to waves resonating 
with our image. In particular the first peak, with the highest 
power, is associated to a wavenumber of fo = 1.59 Gyr, 
and is denoted by the vertical black line in the spectrum. 
The corresponding wavelength of this peak is 0.62 Gyr~^ 
~ 27r/10 Gyr~^. This is in fact the separation we would have 
estimated using Eq. (I1U|I and corresponds to an estimated 
time since accretion of face = 27r/o ~ 10 Gyr . The rest of 
the peaks in the spectrum are associated to the harmonics 
of this wavenumber, i.e., 2/o, 3/o, etc. 

The situation is similar when we consider a cut around 
kr — 0. In this case the scales present in the power spec- 
trum are associated to separations in the direction. We 
find the peak with the highest amplitude to be at the same 
wavenumber fo and it is possible to distinguish some of its 
harmonics. However, due to the smaller number of streams 
we have in the direction, the peaks are less clearly defined 
than for P(kr). 

The Fourier method developed in this section may be 
considered to be very complementary to that introduced by 
MB08. Even if the gravitational potential is only approx- 
imately known, streams in frequency space will still dis- 
tributed on a lattice (see below and Sec. 4.3 of MB08). This 
implies that a Fourier analysis technique will (always) pro- 
vide an estimate of the time since accretion. Instead, in the 
presence of a very clear set of streams MBOS's method can 
be used not only to provide a more accurate estimate of 
the time since accretion but also to narrow down the true 
underlying potential. 



4 FREQUENCY SPACE FOR TIME 
DEPENDENT POTENTIALS 

In a hierarchical universe, the gravitational potentials asso- 
ciated to galaxies are expected to have changed over time. 
Under these circumstances, quantities such as the energy, 
apocentre and pericentre and also the frequencies of indi- 
vidual orbits will no longer be conserved. It is therefore of 
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Figure 6. Top panel: Fourier transform of the distribution of 
streams shown in the top left panel of Fig. [5] The black lines 
denote the regions selected in the Fourier image to compute the 
1-d power spectra, i.e. around fc^ = (middle panel) and fc^ = 
(lower panel). The black vertical lines in these panels show the 
location of the peak with the highest power in the spectrum. 
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Figure 7. Time evolution of the radial action of two different or- 
bits in a very slow (top panel) and a very fast (bottom panel) time 
varying Plummer potentials. The black and light-grey curves cor- 
respond to an outer and an inner orbit respectively. The observed 
oscillatory behaviour is related to the orbital motion. 



special interest to study how a time dependent potential 
afTects the structure of debris in frequency space. 

To this end, we consider a simple model of a satellite ac- 
creted onto a Plummer sphere whose mass varies according 
to 



M{t) = Mo exp{t/ts. 



(17) 



where we explore several values for the time scale, tscaio- 
The growth in mass will lead to a shrinking of the possible 
orbits, and be manifested in an increase in the frequencies 
and acceleration experienced by the particles with time. 



4.1 Time evolution of the frequencies and the 
rate of change of the potential 

To derive the evolution of the frequencies of motion for an 
orbit in a time-dependent mass distribution such as that 
given by Eq. (|17|l . we proceed as follows. At each step of this 
orbit's time- integration, we take the instantaneous position 
and velocity as initial conditions to compute an orbit in the 
instantaneous (frozen) associated potential. Then, for this 
orbit we derive an instantaneous set of frequencies. 

We will consider the evolution of two different orbits: 
an initially outer orbit with an apocentre of rapo ~ 64 kpc 
and pericentre of rper ~ 19 kpc, and a second inner orbit 
with rapo ~ 25 kpc and rper ~ 5 kpc. The initial values of 
the frequencies of the outer orbit are fir « 8 Gyr~^ and 
Sl^ ^ 6 Gyr~^, whereas for the inner one they 
Gyr~^ and fl^ ~ 13 Gyr~^. Both orbits are launched from 
their apocentre. 

Let us first consider a slowly varying potential, where 
iscaie = 30 Gyr. We follow the evolution of the orbit over 
14 Gyr. The top panel in Fig. [7] shows the time evolution of 
the radial action Jr for this case, where black and light-grey 
lines in this Figure correspond to the outer and the inner 
orbit respectively. This Figure shows that the instantaneous 



Figure 8. Time evolution of the radial frequency for two dif- 
ferent orbits in a very slow (top panel) and a very fast (bottom 
panel) time varying Plummer potentials. The black and light-grey 
curves correspond to an outer and an inner orbit respectively. The 
dashed lines represent the time evolution of the normalised ener- 
gies of each orbit E/Eg, while the dotted-dashed line shows the 
time evolution of the mass of the host M/Mq. 



action Jr displays an oscillatory behaviouiQ, with a period 
that corresponds to that of the orbit. 

This is not surprising since under adiabatic varia- 
tions of the potential it is the time average of the actions 
that remains constant (see [Goldstein. Poole fc Safkol 1200 ll : 
Wells & Siklos 2006!). In general, for a given orbit the am- 
plitude of the oscillation depends on the time scale on which 
the potential grows. By increasing fscaie the evolution is more 
adiabatic and the amplitude of the oscillation is decreased. 

In order for an orbit to be in a condition of adiabatic 
invariance it is necessary that tarh ^ iscaie, where torb is the 
orbital period. Since inner orbits have shorter orbital peri- 
ods, we expect those to behave "more adiabatically" com- 
pared to outer orbits. This is demonstrated in the top panel 
of Fig. [T] where the amplitude of the oscillation of Jr is 
smaller for the inner orbit. 

For the outer orbit in the rapidly growing potential, 
iscaic = 3 Gyr, the amplitude of oscillation is larger, espe- 
cially at early times when the orbit is farther away from the 
adiabatic regime (since its initial period, Tr — 0.78, is com- 
parable to tscaic)- As the mass of the system grows in time, 
which leads to the deepening of the potential well, orbits are 
continuously driven towards deeper regions of the well, im- 
plying that the corresponding orbital periods decrease with 
time. Therefore what we observe as the damping of the am- 
plitude of oscillation is a transition from a non-adiabatic 
towards an adiabatic regime. 

In Figure |8] we focus on the evolution of the radial fre- 
quency rir (the angular frequency Sl,^ depicts a qualitatively 
similar behaviour). The top and bottom panels of this Fig- 



^ Note that due to the spherical nature of the potential the other 
two actions, associated with the angular momentum, remain con- 
stant in time. 
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ure correspond, respectively, to the slowly and to the rapidly 
varying potentials. Note that in contrast to the actions, the 
frequencies do change with time, and that their evolution 
appears to follow closely that of the potential. For compar- 
ison, we plot in Fig. |8]the time evolution of the normalised 
energy of each orbit E/Eq (dashed lines), and of the mass 
of the system M/Mq (dotted-dashed line). 

The exact relation between the evolution of the frequen- 
cies and that of the potential can not be derived analytically 
in the general case. However, some insights may be obtained 
from two simple cases: the homogeneous sphere and the Ke- 
pler potentials. 

For the homogeneous sphere the radial period Tr oc p~ 2 
or oc M"5 ijBinnev fc Tremaind [20081 ) . therefore Q.r cx 
M2. Hence, for M{i) cx exp(t/tscaie) we obtain 



Q.r{t) OC exp(t/2tscalc). 



(18) 



Therefore, the frequencies evolve in time with the same func- 
tional form as the potential, but on a time scale that is twice 
as long. Consequently, their evolution is slower. 

In a Kepler potential the radial orbital frequency of 
a particle can be expres sed in terms of its actions as 



ijBinnev fc Tremaindliooa . Appendix E) 



(Jr+LY 



and therefore 



Q,rit) OC exp(2t/f scale) 



(19) 



(20) 



Again, in this example we find that frequencies are 
evolving in a similar fashion as the potential. Therefore a 
very external orbit in the Plummer sphere should initially 
evolve at a rate similar to that given by Eq. (|20|l (as for the 
Kepler potential). On the other hand, inner orbits should 
evolve at a rate similar to that of the homogeneous sphere. 
Thus, we can already see that orbits in a given potential will 
evolve at different rates, depending on their initial location 
in phase-space. 



4.2 Structure in frequency space 

In this section we address how the debris of a satellite is dis- 
tributed in frequency space in the time-dependent Plummer 
potential discussed above. 

We ran simulations with the same satellite and orbital 
initial conditions as described in Section [3.2.21 We followed 
its evolution in the Plummer potential with two different 
timescales, namely tscaic = 3 Gyr and 12 Gyr. As before, we 
do not include the self- gravity of the satellite. 

Fig-Elshows the distribution of particles inside a sphere 
of 4 kpc radius located at « 10 kpc from the centre of the 
system after 10 Gyr of evolution. The top panel presents 
the results for tscaic = 12 Gyr while the bottom panel cor- 
responds to tscaie — 3 Gyr. As we can see from this Figure, 
even in a time dependent system, satellite debris at a given 
spatial location is distributed in a regular pattern of patches 
in the space of orbital frequencies. 

In the simulation where the time evolution is slow (i.e., 
iscaic = 12 Gyr) some of the structures may be decom- 
posed into two clumps. As we explained in Section 13.2.21 
such structures appear when neither the apocentre nor the 
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Figure 9. Distribution of particles in frequency space located 
inside a sphere of 4 kpc radius at 10 kpc from the centre of 
the time-dependent Plummer sphere after 10 Gyr of evolution. 
The timescales for the mass growth of the Plummer sphere are 
iscale = 12 Gyr (top panel) and tgcalo = 3 Gyr (bottom panel). 
Note that structure in frequency space is not destroyed by the 
time varying nature of the potential. 



pericentre of the orbits fall inside the sphere under consid- 
eration. In contrast, when the potential evolves quickly the 
particles in each stream are all distributed in compact sin- 
gle structures. This is because in this case the orbits have 
shrunk so much after 10 Gyr that all streams present in this 
volume have their apocentres inside the chosen sphere. 

Comparison to the case in which the potential is fixed 
(Fig. O shows that streams in the direction are not ex- 
actly distributed along a line of fi^ = est but rather on 
a line with a given curvature. This is due to the fact that 
orbits with different initial frequencies evolve at different 
rates (as discussed in the previous section). The amount of 
curvature observed in this space (e.g. Figure [9} depends on 
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Figure 10. Schematic of the time evolution of the frequencies of 
two particles and their dependence on location in frequency space. 
At initial time, t = to, both particles have the same initial Sl^ but 
different Qr- The frequency of the particle closer to the homoge- 
neous sphere regime evolves at a slower rate than that which is 
closer to the point mass distribution. This results in a lattice that 
has curvature in frequency space for those particles which satisfy 
J^'jsnd Anr(t)(it = 2-rTN, similar to that seen in Figure[9] 

the rate of evolution of the potential. This is schematically 
depicted in Figure [TO] Particles located close to the homo- 
geneous sphere regime, will evolve at rates ~ (2tscaic)~^, 
while those with orbits closer to the point mass distribu- 
tion regime will evolve at rates more similar to ~ 2/fscaic, 
i.e., four times faster. Note that since the particles in the 
Plummer sphere occupy a smaller region in frequency space 
than that delimited by these two regimes, the difference in 
their rate of evolution will typically be smaller than the es- 
timate just provided. If by the final time tend, the particles 
satisfy the condition J^*""'* AQr {t)dt = 2nN their associated 
streams will be located along a curve roughly as depicted in 
this Figure. 

Note that the extent of the distribution of particles in 
frequency space has also changed, but since this is also de- 
pendent on the initial internal properties of the satellite, 
it cannot be used to directly measure the evolution of the 
potential. 

We may conclude from these experiments that the time 
evolution modelled here for the host potential does not de- 
stroy the coherence of satellite debris in frequency space. 
Furthermore, since the frequencies evolve at different rates 
even for particles with same origin, it may be possible to 
recover the evolution of the host potential by measuring 
the curvature of the lines over which the streams are dis- 
tributed in frequency space. Note that this implies that the 
host's evolution has left its imprints in the present day or- 
bits of accreted stars. T his is in contrast to the resu lts by 
iPefiarrubia et al.l (|2006l ) (see also lWarnick et al]|2008l ). who 
find that present-day observables only constrain the present 
day mass distribution of the host, independently of its past 
evolution. 
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Figure 11. Top panel: Fourier transform of the distribution of 
streams shown in the top panel of Fig. [9] As in Fig.|6]the middle 
and central panels depict the 1-d power spectra along the kr and 
directions, respectively. The location of the first peak (denoted 
by the vertical lines in these panels) can be used to estimate the 
accretion time of the satellite. 
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4.3 Estimating the time of accretion 

The top panel in Fig. [11] shows the resuhs of applying a two 
dimensional Fourier transform to an image created from the 
distribution of streams shown in the top panel of Fig. [51 
This corresponds to the satellite accreted 10 Gyr ago onto a 
host whose mass increases on tgcaic = 12 Gyr. As before we 
compute the power spectrum along the and directions. 
The middle panel shows P{kr), and the black line here de- 
notes the location of the peak with the highest amplitude, 
having a wavenumber of /o = 1.46 Gyr. From this spectrum 
we would estimate the satellite was accreted 9.2 Gyr ago. 
When P(k^) is considered, the highest amplitude peak is 
located at a wavenumber /o = 1.36 Gyr, corresponding to 
an estimated time of accretion of 8.6 Gyr. Therefore, in this 
example the analysis of the power spectrum suggests differ- 
ent times of accretion depending on whether il,- or fi^ are 
considered. Both values are reasonably close to the actual 
accretion time (10 Gyr ago), but the direction associated to 
Q.r appears to provide a better estimate. 

Fig. [12] shows the same result but now for the satellite 
accreted onto a host potential that evolves on tscaio = 3 Gyr. 
P{kr) peaks at a wavenumber /o = 1.4 Gyr, as shown by 
the black line in the middle panel of this Figure while P{k^) 
peaks at a wavenumber /o = 1.21 Gyr. The corresponding 
estimated times of accretion tacc are 8.8 and 7.6 Gyr ago, 
respectively. Note that even though the potential evolves on 
a very short timescale, the estimated accretion times only 
differ by ~ 15 — 25% from the true value. 

We may estimate the relation between the true time 
since accretion tacc and tacc for a general time dependent po- 
tential as follows. Let us consider two particles that at final 
time are found in adjacent streams. As discussed before, the 
condition for this to happen is A6i(tacc) = Jo'^"" AD,i{t)dt = 
2n, for i = r, 0. This implies 



Jo 



(21) 



We may use the se cond mean value theorem for integra- 
tion (|Courantlll99i ). i.e. G{t)-f{t)dt = G(a) <fi{t)dt + 
G{b) J^ip{t)dt. If we take ip{t) = 1 and G{t) = Ar2i(t), then 

/ Ailiit) dt = AQi{0)tg + Ani(facc)(iacc - tg), 

Jo 

where tg £ (0, face). Therefore in Eq. ([gTJ 



ta 



face - tg(l - Af7,(0)/Af7,(facc)), 



(22) 



which shows that our estimate is always a lower limit to 
the actual time since accretion. Its accuracy depends on the 
particular form of the gravitational potential and on its evo- 
lution (through the quantities Ani(tacc) and t^). 

We have also performed an experiment in which an ax- 
isymmetric thin disk was grown adiabatically inside an ini- 
tially spherical (Plummer) mass distribution. Therefore, in 
this case, both the shape and the depth of the gravitational 
potential have changed with time. We find that also under 
such conditions, we are able to estimate the time of accretion 
with similar accuracy. 
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Figure 12. As in Fig. 1111 now for the distribution of streams 
shown in the bottom panel of Fig. |9] 
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5 FULL N-BODY CASE 

The previous analysis has shown that frequency space is 
particularly well suited for identifying streams from merg- 
ers, at least for the idealised potentials considered thus far. 
We now explore a more realistic case, namely that of a mi- 
nor merger with a live hos t. We use one of the simulations of 
IVillalobos fc Helmil l|2008l ) (hereafter VH08), that model the 
formation of a thick disk via the merger of a relatively mas- 
sive satellite onto a pre-existing thin disk. Two important 
physical processes take place during the merger which may 
affect the structure of the debris in frequency space. First, 
the satellite suffers significant dynamical friction. Secondly, 
the host system responds strongly to the perturber, result- 
ing in a disk that has been significantly tilted and heated. 
This implies that, even though the total mass of the sys- 
tem is conserved, its distribution evolves in time. The VH08 
simulation that we consider here corresponds to a two com- 
ponent satellite (stars -I- dark matter) launched on a 30 ° 
inclination orbit with respect to the disk from a distance of 
approximately 84 kpc (~ 50 disk scale lengths). The host 
consists of a (live) dark halo and a (live) thin disk. The to- 
tal (and the stellar) mass ratio between the satellite and the 
host is 20%. After 4.5 Gyr of evolution, the satellite is fully 
disrupted and has deposited debris in the disk of the host. 
The host disk is significantly thickened by this process and 
shows characteristics typical of thick disks. 

5.1 Computing the frequencies 

The remnant system modelled by VH08 does not have a 
separable gravitational potential which prevent us from de- 
riving explicitly the action-angles variables using Eq. (|2])- 
([3]). Nevertheless several methods have been developed to 
obtain the frequencies of orbits in general potentials. Here 
we focus on th e spectral analysis approach, introduced by 
iBinnev fc Sperge l (1982), which relies on the numerical in- 
tegration of the equations of motion in a given potential. 
The basic idea is to perform a Fourier transform of the time 
series x(t) (i.e., the orbit), and then to derive the frequencies 
of motion. 

To obtain a reliable Fourier spectrum the orbits must 
be integrated for a very large number of periods. However, 
in the problem under consideration the orbits are not sta- 
tionary, i.e. the frequencies are time dependent and so this 
approach cannot be applied directly using the orbits from 
the N-body simulation. This is why we have used the final 
output of the simulation to define a set of initial conditions 
which we integrate in a fixed (static) potential. This poten- 
tial should ressemble that of our simulation, and for simplic- 
ity we have taken a t wo component model, con sisting in a 
Miyamoto-Nagai disk (iMivamoto fc Nagaiiri975h 



^R^ + {a + v/i2-rF)2 

and NEW dark matter halo (|Navarro. Frenk fc Whit^llQQel ) 

4'haio = -<&o-logfl + -V (24) 

We have set the parameters of the model to Afdisk = 2.5 x 
lOi°M0, a = 1.8 kpc and b = 0.6 kpc; <l>o = 1.05 x 10^ 
km^ s~^ and Vs — 14.1 kpc. These choices lead to a potential 



energy distribution that differs by at most 5 per cent from 
the true potential up to a radius of 50 kpc on the z = Q plane. 
This radius is large enough to enclose the apocentres of 99% 
of the particles presently found in a volume resembling the 
Solar neighbourhood. We thus integrate the orbits of the 
stellar particles located in such a volume for approximately 
100 orbital periods and compute their frequenc ies using the 
code developed bv lCarpintero fc Aguilaii (|l998r ). 

Note that the results are not strongly dependent on our 
particular choice of th e potential, since substructure in inte- 
grals of motion space (|Helmi fc de Zeeuwll200ol : iHelmi et al.l 
2006) and in frequency space (MB08) are both robust to 
small differences in the mass distribution. This is because 
we always focus on small volumes in space and hence small 
changes in the potential essentially act as a zero point off- 
set, affecting all particles present in this volume in the same 
way. 

5.2 Results 

In Fig. [13] we present the distribution of stellar particles in- 
side a sphere of 2 kpc radius located at 8 kpc from the cen- 
tre of the remnant galaxy in some commonly used spaces to 
identify substructure. Here the red particles represent those 
originally present in the disk, while those from the satellite 
are colour-coded in black. The total number of stellar parti- 
cles inside this sphere is approximately 27000, out of which 
23000 belong to the satellite and 4000 to the hosfl 

The top panel of Fig. [TH] shows the distribution of 
particles in the Q.^ vs. Q.r space. Note the large amount 
of substructure present. As expected, satellite particles are 
distributed in multiple lumps associated with the different 
streams crossing this volume. These lumps are better de- 
fined for orbits which spend most of their time far from the 
centre, and so are found at low frequencies and have long 
mixing timescales. 

The upper envelope defined by the smallest Q.r at a 
given Q.^ corresponds to particles with low eccentricity (i.e. 
on more circular orbits) as discussed in Section 13.2.11 In 
our simulation this region is preferentially populated by disk 
particles. 

In general, orbits in axisymmetric potentials have three 
independent frequencies: the previously discussed Q.r and Q.^ 
and a third frequency associated with the vertical motion, 
Slz. However, particles in the thick disk typically have very 
short periods in the vertical direction. Hence, the number 
of streams expected in the fiz direction is very large and 
therefore it is hard to resolve in our simulation, and not 
discussed further here. 

In the bottom left panel of Fig. [13] we can observe the 
distribution of particles in the space defined by orbital apoc- 
entre and pericentre. The stripy features along lines of differ- 
ent eccentricities are the counterparts of the lumps located 
on diagonal lines in frequency space. The middle panel shows 
the E — Lz space, where extended structures with slightly 
different orientations may be seen. The panel on the right 

■2 This is due to higher number of particles used to simulate the 
satellite to properly resolve its structure in phase-space. There- 
fore, in practice only 1 in 50 satellite particles should be consid- 
ered. 
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Figure 13. Distribution of stellar particles inside a sphere of 4 kpc radius located at 8 kpc from the centre of the remnant of a 5:1 
merger between a satellite and a disk galaxy, after 4.5 Gyr of evolution. In all panels, the black dots represent particles from the satellite 
and red dots particles from the disk. The top panel shows the distribution of particles in frequency space, whereas the bottom panels 



correspond to the apocentre vs. pericentre space (left), the E — Lz (middle) and the 
are much better defined in frequency space than in these other projections. 



(right) projections. Notice that streams 



shows a projection of velocity space, which clearly exhibits 
less well-defined structures in comparison to the other spaces 
discussed. In general, we may conclude from this Figure that 
streams are much better defined and easier to interpret in 
frequency space than in any of the other projections consid- 
ered. 



5.3 Estimating the time of accretion 

In comparison with the idealised cases discussed in Section 
I3.2.2l and l4.2l the distribution of particles in frequency space 
for the VII08 simulation is less regular. Nevertheless, we 
would like to explore if it is still possible to obtain a good 
estimation of the time of accretion of the satellite, using the 
Fourier analysis described in previous sections. 

The top panel of Fig. [14] plots the image in Fourier 
space of the frequency plane Q,^ vs. fir discussed above. This 
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image was obtained by making a grid in frequency space 
and counting the number of particles in each grid element. 
As before the black lines in the top panel of Fig. 1141 denote 
the cuts around kr = and = made to compute the 
1-dimensional power spectra. 

To estimate the uncertainties in the power spectrum, 
we create a second image obtained by assuming that the 
particles are spread throughout the permitted regions of 
frequency space according to a Poisson distribution. That 
is, we assign to each grid element A'j, particles, where this 
number is drawn from a Poisson distribution with mean 
(A'") = (Adisk + Asat)/wgrid, i.e. equal to average number 
of stellar particles per grid element. We then compute the 
Fourier transform of this distribution and the corresponding 
1-dimensional power spectra from the image. This procedure 
is repeated 1000 times, to yield an average power spectrum 
in each direction associated to a random distribution of par- 
ticles in frequency space. 

The results are shown in the bottom and middle panel of 
Fig. 1141 The black solid line in each panel denotes the power 
spectrum obtained from the actual distribution of particles 
in frequency space. The dashed line corresponds to the av- 
erage power spectrum for our random distributions where 
the vertical error bars represent the rms dispersion around 
this average spectrum. As we can see, a statistically signifi- 
cant peak can be identified in the observed power spectrum 
along kr (and also for fc^). The wavenumber of this peak is 
/o ~ 0.43 Gyr for both the radial and angular directions. 
This corresponds to an estimated accretion time of 2.7 Gyr. 
Recall that this simulation is evolved for 4.5 Gyr. However, 
the satellite only starts to lose a significant amount of stars 
~ 1 Gyr after infall (see Figure 3 of VH08). Therefore, we 
may conclude that the time since disruption of a satellite 
can be reasonably estimated even in a simulation with a live 
host. 
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Figure 14. Top panel: Fourier transform of the distribution of 
particles in frequency space shown in the top panel of Fig. 1131 
As in Fig. |6] the middle and central panels depict the 1-d power 
spectra around fc^ ~ and kr ~ 0, respectively. The vertical 
lines in these panels show the location of the peak associated 
to the average separation between streams in frequency space, 
and can be used to estimate the accretion epoch. The dashed 
curve corresponds to the average power spectrum obtained from 
a random distribution of particles in frequency space. The error 
bars represent the rms dispersion from this average spectrum. 



6 DISCUSSION AND CONCLUSIONS 

We have confirmed and extended previous work by MB08 
showing that orbital frequencies constitute a very suitable 
space to identify debris from past accretion events. In this 
space, particles in a given volume of physical space are found 
to populate well-defined lumps, each of them associated with 
a different stream. 

For time-independent potentials, these streams are dis- 
tributed in a regular pattern along lines of constant fre- 
quency. As time goes by, the number of streams at a given 
physical location increases while the separation between ad- 
jacent streams decreases. We have shown here that this char- 
acteristic separation or scale can be used to estimate the 
time of accretion of the object through a Fourier analysis. 

We have also addressed how the time evolution of the 
host affects substructure in frequency space. As an example, 
we have considered the case in which the mass of the host is 
increased exponentially in time on two different timescales. 
We find that in contrast with the actions which are adiabatic 
invariants for slowly varying potentials, the frequencies al- 
ways evolve in time, closely following the rate of change of 
the host potential. This evolution however does not destroy 
the dumpiness present in frequency space. Streams still look 
lumpy and are regularly distributed in this space, even if the 
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host potential varies on a very short timescale (when even 
the actions are no longer invariant). Interestingly in this 
case, streams are not exactly distributed along straight lines 
of constant frequency but rather on lines whose curvature 
depends on the timescale of growth of the potential. This im- 
plies that, contrary to prev ious claims fe.g. lPenarrubia et al.l 
I2OO6I : IWarnick et al.ll2008l ). the final distribution of streams 
does retain information on the evolution in time of the host. 

Finally, we have analysed a full N-body simulation of 
the accretion of a satellite galaxy onto a disk galaxy. Due to 
the inclusion of other physical processes (such as self-gravity, 
dynamical friction and the variation in time of the distribu- 
tion of mass in the system) , the distribution of satellite par- 
ticles in frequency space looks somewhat less regular than 
in the previously discussed idealised cases. Nevertheless, the 
space of frequencies is still rich in substructure associated 
to streams. Furthermore, even in this case Fourier analysis 
techniques can be used to estimate the time of accretion 
reasonably well, although in general, only a lower limit is 
obtained. 

For all these examples we have compared the final distri- 
bution of streams in the most commonly used spaces to iden- 
tify satellite debris. These comparisons have shown that fre- 
quency space contains information that is either not present 
or is not simply obtained in other spaces. In general streams 
are most sharply defined in frequency space. 

One important simplification in our analysis is to con- 
sider the accretion of a single satellite galaxy. In reality, the 
process of the formation of a galaxy like the Milky Way 
will have involved many me rgers (in the range of 5 — 40, 
see iDe Lucia fc Helmill2008l ). Therefore it is possible, and 
even likely, that their debris will overlap in frequency space. 
However, as we have shown, the separation between adja- 
cent streams in frequency space is predicted to be different 
for mergers that took place at different epochs. Moreover, 
the location of their debris in frequency space at the present 
time is dependent upon the initial orbital conditions. Con- 
sequently we expect the overlap to be incomplete and hence 
the identification of their remnants to be feasible. We are 
currently analysing the full cosmological hig h resolution N- 
body simulations of the Aquarius Project (jSpringel et al.l 
|2008| ) in order to address this issue and to understand what 
to expect in the context of the hierarchical paradigm of 
structure formation. 
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